%%
Ts = [12,15,20,23,26,30,35,40,50,60,80,100,130,160,190,220,250,275,296];
dataall=zeros(19,127);
files = ['HV';'VH'];
pol = ['HV';'VH'];
fields = 0:0.5:7;

bs=[510,514,518];

%eqn='a1/(1+(x-b1)^2/c1^2)+a2/(1+(x-b2)^2/c2^2)+a3/(1+(x-b3)^2/c3^2)+a4/(1+(x-b4)^2/c4^2)+a5/(1+(x-b5)^2/c5^2)+e*x+f';
%sp=[10,30,50,30,10, bs(1), bs(2), bs(3), bs(4), bs(5),.4,.4,.4,.4,.4,0,0]
eqn='a1/(1+(x-b1)^2/c1^2)+a2/(1+(x-b2)^2/c1^2)+a3/(1+(x-b3)^2/c1^2)+e*x+f';
sp=[10,30,200, bs(1), bs(2), bs(3),.4,0,0]
eqn2='a2/(1+(x-b2)^2/c1^2)+a3/(1+(x-b3)^2/c1^2)+e*x+f';
sp2=[30,200, bs(2)-1, bs(3)-1,.4,0,0]
eqn3='a3/(1+(x-b3)^2/c1^2)+e*x+f';
sp3=[100, bs(3)-3,2,0,0]
%eqn='a1*exp(-(x-b1)^2/c1^2)+a2*exp(-(x-b2)^2/c2^2)+a3*exp(-(x-b3)^2/c3^2)+a4*exp(-(x-b4)^2/c4^2)+a5*exp(-(x-b5)^2/c5^2)+a6*exp(-(x-b6)^2/c6^2)+e*x+f';
%sp=[10,30,50,30,10, bs(1), bs(2), bs(3), bs(4), bs(5),.4,.4,.4,.4,.4,0,0]
%eqn='a1*exp(-(x-b1)^2/c1^2)+a2*exp(-(x-b2)^2/c1^2)+a3*exp(-(x-b3)^2/c1^2)+a4*exp(-(x-b4)^2/c1^2)+a5*exp(-(x-b5)^2/c1^2)+a6*exp(-(x-b6)^2/c1^2)+e*x+f';
%sp=[10,30,50,30,10,10, bs(1), bs(2), bs(3), bs(4), bs(5), bs(6),.4,0,0]
b1s=[];
b2s=[];
b3s=[];
b1serr=[];
b2serr=[];
b3serr=[];
c1s=[];
c1serr=[];



for filenum=2:2
    b4s=[];
    figure(23)
    hold on
    for num=1:19
        filename=[files(filenum,:),num2str(Ts(num)),'.txt'];
        data=readmatrix(filename);
        Es=data(1779:1905,1);
        Is=data(1779:1905,2);
        if num<12.5
            fo = fitoptions('Method','NonlinearLeastSquares',...
            'StartPoint',sp);
            ft = fittype(eqn,'options',fo);
            [curve2,gof2] = fit(Es,Is,ft);
            b1s=[b1s,curve2.b1];
            b2s=[b2s,curve2.b2];
            b3s=[b3s,curve2.b3];
            c1s=[c1s,curve2.c1];
            conf = confint(curve2, 0.95);
            standardErrors = (conf(2,:) - conf(1,:)) / 2;
            b1serr=[b1serr,standardErrors(4)];
            b2serr=[b2serr,standardErrors(5)];
            b3serr=[b3serr,standardErrors(6)];
            c1serr=[c1serr,standardErrors(7)];
        end
        if num>12.5 && num<15.5
            fo = fitoptions('Method','NonlinearLeastSquares',...
            'StartPoint',sp2);
            ft = fittype(eqn2,'options',fo);
            [curve2,gof2] = fit(Es,Is,ft);
            b2s=[b2s,curve2.b2];
            b3s=[b3s,curve2.b3];
            c1s=[c1s,curve2.c1];
            conf = confint(curve2, 0.95);
            standardErrors = (conf(2,:) - conf(1,:)) / 2;
            b2serr=[b2serr,standardErrors(3)];
            b3serr=[b3serr,standardErrors(4)];
            c1serr=[c1serr,standardErrors(5)];
        end
        if num>15.5
            fo = fitoptions('Method','NonlinearLeastSquares',...
            'StartPoint',sp3);
            ft = fittype(eqn3,'options',fo);
            [curve2,gof2] = fit(Es,Is,ft);
            b3s=[b3s,curve2.b3];
            c1s=[c1s,curve2.c1];
            conf = confint(curve2, 0.95);
            standardErrors = (conf(2,:) - conf(1,:)) / 2;
            b3serr=[b3serr,standardErrors(2)];
            c1serr=[c1serr,standardErrors(3)];
        end
        dataall(num,:)=Is-curve2.e.*(Es)-curve2.f;
        plot(data(:,1),data(:,2)+(num-1)*50-curve2.e.*(data(:,1))-curve2.f,'.','Color',[(num)/19 0 0]);
        plot(505:0.02:525,curve2(505:0.02:525)+(num-1)*50-curve2.e.*(505:0.02:525)'-curve2.f,'-','Color',[(num)/19 0 0]);
        %b4s(num)=curve2.c1;
    end
    ylim([-20,1000])
    xlim([505,525])
    %legend('RL','RR','LL','LR')
    title(pol(filenum,:))
    %}
    %%figure
    %imagesc(Es,fields,dataall)
    %colormap('jet')
    %clim([20,80])
    %%[X,Y]=meshgrid(Es,fields);
    %%mesh(X,Y,dataall,'LineWidth',1.5)
    %%clim([-0,50])
    %%xlim([285,310])
    %%colormap('jet')
    %%title(files(filenum,:))
    %%figure(200)
    %%hold on
    %%plot(fields,b4s)
end
eqnfit1='-a*(1+2/(exp((c)/(2*0.69501*x))-1))+c';
spfit1=[3,522];
spfit2=[3,517];
spfit3=[3,513];
fo = fitoptions('Method','NonlinearLeastSquares',...
            'StartPoint',spfit1);
ft = fittype(eqnfit1,'options',fo);
[curve1,gof2] = fit(Ts',b3s',ft);
fo = fitoptions('Method','NonlinearLeastSquares',...
            'StartPoint',spfit1);
ft = fittype(eqnfit1,'options',fo);
[curve2,gof2] = fit(Ts(1:15)',b2s',ft);
fo = fitoptions('Method','NonlinearLeastSquares',...
            'StartPoint',spfit1);
ft = fittype(eqnfit1,'options',fo);
[curve3,gof2] = fit(Ts(1:12)',b1s',ft);
figure(50);hold on;
errorbar(Ts,b3s,b3serr,'o')
plot(0:300,curve1(0:300))
errorbar(Ts(1:15),b2s,b2serr,'o')
plot(0:300,curve2(0:300))
errorbar(Ts(1:12),b1s,b1serr,'o')
plot(0:300,curve3(0:300))
%

figure(51);hold on;
errorbar(Ts,c1s,c1serr,'o')
%% Imagesc plot
Tplot=0:305;
dataplot=zeros(length(Tplot),127);
Tline=[];
for num=1:(length(Ts)-1)
    Tline(num)=(Ts(num)/2+Ts(num+1)/2)-0.1;
end
for idx=1:length(Tplot)
    numT=sum(Tplot(idx)>Tline)+1
    dataplot(idx,:)=dataall(numT,:);
end
figure(52);
imagesc(Es,Tplot,dataplot)
colormap("jet")
clim([0,80])
xlim([500,530])
ylim([0,300])
%%

dataall=zeros(15,127);
files = ['HV';'VH'];
pol = ['LL';'RR'];
fields = 0:0.5:7;

bs=[356.95,361.484,363.423];

%eqn='a1/(1+(x-b1)^2/c1^2)+a2/(1+(x-b2)^2/c2^2)+a3/(1+(x-b3)^2/c3^2)+a4/(1+(x-b4)^2/c4^2)+a5/(1+(x-b5)^2/c5^2)+e*x+f';
%sp=[10,30,50,30,10, bs(1), bs(2), bs(3), bs(4), bs(5),.4,.4,.4,.4,.4,0,0]
eqn='a1/(1+(x-b1)^2/c1^2)+a2/(1+(x-b2)^2/c1^2)+a3/(1+(x-b3)^2/c1^2)+e*x+f';
sp=[10,30,200, bs(1), bs(2), bs(3),.4,0,0]
%eqn='a1*exp(-(x-b1)^2/c1^2)+a2*exp(-(x-b2)^2/c2^2)+a3*exp(-(x-b3)^2/c3^2)+a4*exp(-(x-b4)^2/c4^2)+a5*exp(-(x-b5)^2/c5^2)+a6*exp(-(x-b6)^2/c6^2)+e*x+f';
%sp=[10,30,50,30,10, bs(1), bs(2), bs(3), bs(4), bs(5),.4,.4,.4,.4,.4,0,0]
%eqn='a1*exp(-(x-b1)^2/c1^2)+a2*exp(-(x-b2)^2/c1^2)+a3*exp(-(x-b3)^2/c1^2)+a4*exp(-(x-b4)^2/c1^2)+a5*exp(-(x-b5)^2/c1^2)+a6*exp(-(x-b6)^2/c1^2)+e*x+f';
%sp=[10,30,50,30,10,10, bs(1), bs(2), bs(3), bs(4), bs(5), bs(6),.4,0,0]
fo = fitoptions('Method','NonlinearLeastSquares',...
        'StartPoint',sp);

for filenum=1:2
    b4s=[];
    figure(24)
    hold on
    for num=1:1
        filename=[files(filenum,:),num2str(fields(num)),'T.txt'];
        data=readmatrix(filename);
        Es=data(1779:1905,1);
        Is=data(1779:1905,2);
        ft = fittype(eqn,'options',fo);
        [curve2,gof2] = fit(Es,Is,ft);
        dataall(num,:)=Is-curve2.e.*(Es)-curve2.f;
        plot(data(:,1),data(:,2)+(filenum-1)*0-curve2.e.*(data(:,1))-curve2.f,'.','Color',[(num)/15 0 0]);
        plot(490:0.02:530,curve2(490:0.02:530)+(filenum-1)*0-curve2.e.*(490:0.02:530)'-curve2.f,'-','Color',[(num)/15 0 0]);
        b4s(num)=curve2.c1;
    end
    ylim([-50,700])
    xlim([350,370])
    %legend('RL','RR','LL','LR')
    title(pol(filenum,:))
    %}
    %%figure
    %imagesc(Es,fields,dataall)
    %colormap('jet')
    %clim([20,80])
    %%[X,Y]=meshgrid(Es,fields);
    %%mesh(X,Y,dataall,'LineWidth',1.5)
    %%clim([-0,50])
    %%xlim([285,310])
    %%colormap('jet')
    %%title(files(filenum,:))
    %%figure(200)
    %%hold on
    %%plot(fields,b4s)
end

